In [1]:
from util import *
from glob import glob

/usr/local/lib/python3.10/dist-packages/geopandas/_compat.py:123: UserWarning: The Shapely GEOS version (3.11.1-CAPI-1.17.1) is incompatible with the GEOS version PyGEOS was compiled with (3.10.4-CAPI-1.16.2). Conversions between both will be slow.
  warnings.warn(
/home/nyou045/git/CliffErosion_Projections/util.py:3: UserWarning: Shapely 2.0 is installed, but because PyGEOS is also installed, GeoPandas will still use PyGEOS by default for now. To force to use and test Shapely 2.0, you have to set the environment variable USE_PYGEOS=0. You can do this before starting the Python process, or in your code before importing geopandas:

import os
os.environ['USE_PYGEOS'] = '0'
import geopandas

In a future release, GeoPandas will switch to using Shapely by default. If you are using PyGEOS directly (calling PyGEOS functions on geometries from GeoPandas), this will then stop working and you are encouraged to migrate from PyGEOS to Shapely 2.0 (https://shapely.readthedocs.io/en/latest/migration_pygeos.html).
  import geopandas as gpd  # Geospatial data
In [2]:
df = load_AOIs()
df
Out[2]:
Taranaki AOI SSP 4.5 (p50) SSP 4.5 (p83) SSP 8.5 (p50) SSP 8.5 (p83) Rate SSP 4.5 (p50) Rate SSP 4.5 (p83) Rate SSP 8.5 (p50) Rate SSP 8.5 (p83) match match_score
7 NORTH TongaporutuRiver 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 TongapurutuRiverCliffs 93.750000
11 SOUTH HangatahuaRiver_South 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 HangatahuRiver_South 97.435897
21 SOUTH Rahotu 0.58 0.78 0.84 1.10 0.0058 0.0078 0.0084 0.0110 Rahotu 100.000000
20 SOUTH Pihama 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Pihama 100.000000
19 SOUTH OpunakeBeach 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 OpunakeBeachCliffs 100.000000
18 SOUTH OhaweBeach 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 OhaweBeach 100.000000
17 SOUTH Oeo 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Oeo 100.000000
16 SOUTH Manutahi 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 Manutahi 100.000000
15 SOUTH ManaBay 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 ManaBayCliffs 100.000000
14 SOUTH KaupokonuiBeach 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 KaupokonuiBeach 100.000000
13 SOUTH Kakaramea 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 Kakaramea 100.000000
12 SOUTH Hawera_WaihiBeach 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 Hawera_WaihiBeach 100.000000
0 NORTH Mohakatino_River 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 MohakatinoRiver 100.000000
10 SOUTH CapeEgmont 0.58 0.78 0.84 1.10 0.0058 0.0078 0.0084 0.0110 CapeEgmont 100.000000
9 NORTH Waitara 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Waitara 100.000000
8 NORTH UrenuiRiver_North 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 UrenuiRiverNorthCliffs 100.000000
6 NORTH PariokariwaPoint 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 PariokariwaPointCliffs 100.000000
5 NORTH Onaero 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 OnaeroCliff 100.000000
4 NORTH Oakura_South 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 OakuraSouth 100.000000
3 NORTH Oakura 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Oakura 100.000000
2 NORTH NewPlymouth_Airport 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 NewPlymouthAirport 100.000000
1 NORTH NewPlymouth 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 NewPlymouthCliffs 100.000000
22 SOUTH WainuiBeach 0.57 0.78 0.83 1.09 0.0057 0.0078 0.0083 0.0109 WainuiBeach 100.000000
23 SOUTH WaipipiBeach 0.57 0.78 0.83 1.09 0.0057 0.0078 0.0083 0.0109 WaipipiBeachCliffs 100.000000
In [3]:
site = df.match.sample(1).iloc[0]
site
Out[3]:
'Oakura'
In [4]:
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
if site == "ManaBayCliffs":
  print("Flipping")
  for k, v in transect_metadata.items():
    transect_metadata[k]["Azimuth"] = v["Azimuth"] + 180
In [5]:
results = predict(gdf, transect_metadata)
results
Out[5]:
TransectID BaselineID group Year ocean_point linear_model_point sqrt_model_point BH_model_point Sunamura_model_point
0 8 1 0 2100 POINT (1682577.3607581342 5671120.99636598) POINT (1682903.9877330572 5670743.528687456) POINT (1682904.7560665228 5670742.640760029) POINT (1682880.738653313 5670770.396572942) POINT (1682903.9877330572 5670743.528687456)
1 9 1 0 2100 POINT (1682590.3312759965 5671130.686919325) POINT (1682916.521771773 5670744.2823856585) POINT (1682927.4639775786 5670731.32027307) POINT (1682893.6026720947 5670771.432300012) POINT (1682916.521771773 5670744.2823856585)
2 10 1 0 2100 POINT (1682598.6770314383 5671138.73452017) POINT (1682927.0276047606 5670747.144026966) POINT (1682945.154076735 5670725.5264186) POINT (1682904.1986634596 5670774.369794513) POINT (1682927.0276047606 5670747.144026966)
3 11 1 0 2100 POINT (1682604.6288549996 5671145.294702299) POINT (1682937.4481723008 5670749.638776907) POINT (1682964.3080712678 5670717.70770311) POINT (1682914.576493046 5670776.828651263) POINT (1682937.4481723008 5670749.638776907)
4 12 1 0 2100 POINT (1682607.6075595547 5671154.595712404) POINT (1682942.7813369657 5670757.408985613) POINT (1682972.087082433 5670722.681170065) POINT (1682919.8669761904 5670784.56289969) POINT (1682942.7813369657 5670757.408985613)
... ... ... ... ... ... ... ... ... ...
97 187 4 7 2100 POINT (1685286.8791120881 5673060.735761237) POINT (1685696.9093590495 5672775.544998921) POINT (1685703.3645548837 5672771.055178392) POINT (1685667.740752512 5672795.832812311) POINT (1685696.9093590495 5672775.544998921)
98 188 4 7 2100 POINT (1685286.8521959984 5673068.408532419) POINT (1685700.9186097626 5672784.572777532) POINT (1685707.2977494132 5672780.199981628) POINT (1685671.612585248 5672804.661578058) POINT (1685700.9186097626 5672784.572777532)
99 189 4 7 2100 POINT (1685287.0361195225 5673075.640930105) POINT (1685706.5802912419 5672792.548262031) POINT (1685720.9431551918 5672782.856740254) POINT (1685677.127796816 5672812.421700373) POINT (1685706.5802912419 5672792.548262031)
100 190 4 7 2100 POINT (1685289.1618428899 5673080.93252381) POINT (1685743.416889592 5672780.039583565) POINT (1685834.9939295463 5672719.380081791) POINT (1685713.7955143917 5672799.660417153) POINT (1685743.416889592 5672780.039583565)
101 191 4 7 2100 POINT (1685294.8243254684 5673083.481260104) POINT (1685719.8676822786 5672807.560334522) POINT (1685733.991497258 5672798.391726098) POINT (1685690.0660794459 5672826.906327597) POINT (1685719.8676822786 5672807.560334522)

102 rows × 9 columns

In [6]:
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
poly = prediction_results_to_polygon(results)
In [7]:
poly.explore(tiles="Esri.WorldImagery")
Out[7]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [34]:
site = df.match.sample(1).iloc[0]
site = "CapeEgmont"
print(site)
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
m = gpd.read_file(f"Shapefiles/{site}_TransectLines.shp").explore("TransectID", tiles="Esri.WorldImagery")
results = predict(gdf, transect_metadata)
results["geometry"] = results["linear_model_point"]
pd.concat([gdf, results])[["Year", "geometry"]].explore("Year", m=m)
gpd.GeoSeries(results.ocean_point, crs=2193).explore(m=m, color="blue")
#gpd.GeoSeries(poly, crs=2193).explore(color="pink", m=m)
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
polygon = prediction_results_to_polygon(results)
polygon.explore(m=m)
m

CapeEgmont
Out[34]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [9]:
run_all_parallel()

  0%|          | 0/24 [00:00<?, ?it/s]
Flipping
In [10]:
samples = pd.concat(gpd.read_file(f) for f in glob("Projected_Shoreline_Polygons/*_linear_Rate SSP 4.5 (p50).shp"))
samples.explore(tiles="Esri.WorldImagery")
Out[10]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [11]:
# Compare different SLR rates
site = df.match.sample(1).iloc[0]
print(site)
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
low_proj_slr = predict(gdf, transect_metadata, Proj_SLR=0.007)
high_proj_slr = predict(gdf, transect_metadata, Proj_SLR=0.015)
m = gpd.GeoSeries(low_proj_slr.sqrt_model_point, crs=2193).explore(color="blue", tiles="Esri.WorldImagery")
gpd.GeoSeries(high_proj_slr.sqrt_model_point, crs=2193).explore(color="green", m=m)

WaipipiBeachCliffs
Out[11]:
Make this Notebook Trusted to load map: File -> Trust Notebook